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Abstract 

We consider a pair of stochastic integrate and fire neurons receiving correlated stochastic inputs. The evolution of 
this system can be described by the corresponding Fokker-Planck equation with non-trivial boundary conditions re- 
sulting from the refractory period and firing threshold. We propose a finite volume method that is orders of magnitude 
faster than the Monte Carlo methods traditionally used to model such systems. The resulting numerical approxima- 
tions are proved to be accurate, nonnegative and integrate to 1. We also approximate the transient evolution of the 
system using an Ornstein-Uhlenbeck process, and use the result to examine the properties of the joint output of cell 
pairs. The results suggests that the joint output of a cell pair is most sensitive to changes in input variance, and less 
sensitive to changes in input mean and correlation. 



1 Introduction 



The integrate and fire (IF) model is used widely in mathematical biology ( Burkitt 2006 , Keener and Sneyd 2008| l. It is 
simple, yet versatile, and provides a good approximation of the response of an excitable cell in a variety of situations. 
A stochastic version of the IF model can describe the behavior of large populations of cells through the evolution of 
the corresponding probability density ( Kni ght] |1972||Nykamp and Tranchina||2000||Rolls et aL]|2008| l. It can also be 



used to study the response of a single cell subject to a large number of small, statistically independent inputs ( |Lindner 
[200T] |Renart etaLl [2003) 

Collections of excitable cells frequently do not behave independently. The joint response of populations of electri- 
cally active cells is of interest in a number of areas in biology: Pancreatic j3 -cells have to synchronize their response to 
secrete insulin (Meda et al. 1984 Sherman and Rinzel||1991[ l, and the coordinated activity of cardiac cells is essential 
for their function (Keener and Sneyd 2008). Our study is motivated primarily by cells in neural populations. Such 
cells typically fire action potentials in response to synaptic inputs from other cells. The standard stochastic IF model 



can capture the response of a cell when such inputs are independent (Renart et al. 2003 i. However dependencies 
between these inputs cannot always be ruled out. Such dependencies can affect the output statistics of a neuronal 
population, and significantly impact the amount of information carried in the population response (Salinas a nd Se-| 
jnowski 2000[ Sompolinsky et al. 200 1 ji . Even weak correlati ons between individual cells can significantly impact 



the ensemble activity of a population (Shadlen and Newsome 1998[ Rosenbaum et al. 2010| l. Here we describe a 



model that can be used to examine the behavior of two cell populations (or cell pairs) receiving correlated inputs. 

We first develop a Fokker-Planck equation that describes the evolution of the probability density for a pair of cells 
receiving correlated inputs. The response of cell pairs receiving correlated inputs has been studied previously using 
linear response theory (Ostoj ic~et al.||20 09; de la Roc ha et aL] |2"007), and numerical simulations (Galan et al, 2007) in 
related models. However, the boundary conditions in the presence of a refractory period are nontrivial, and can impact 
the behavior of the system. We therefore present the model in some detail. 

Next we describe a finite volume method that can be used to study solve the Fokker-Planck equation numerically 
for the probability density. Previously, we proposed a fast and accurate finite volume method for modeling a general IF 
neuron driven by a stochastic input ( Marpeau et al. , 2009). This method was significantly faster than Monte Carlo (MC) 
simulations and we proved several stability properties of the algorithm. Here we extend this method to two neurons 
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with correlated inputs. While the dynamics of interacting populations has been examined previously (Nyka mp and| 
|Tranchina] |2000| |Harrison| |et al 2 005 ), we are not aware of a numerical treatment of the Fokker-Planck equation 
corresponding to stochastic IF neurons driven by correlated noise. 

Finally, we develop a simple analytical approximation in terms of a related Ornstein-Uhlenbeck process that cap- 
tures the response of a cell pair to study the behavior of cells, or cell populations receiving correlated inputs. This 
approach provides an alternative to the linear response techniques commonly in use fLindner| |2001| |Ostojic et al.| 
|2D09l >. 

We use this approximation to examine the response of a single cell and a cell pair to changes in the input parameters. 
The variance (noisiness), mean and synchrony between the inputs are separate channels along which information can 
be communicated to postsynaptic cells. We find that the spiking statistics of a single cell and the cell pair are most 
sensitive to changes in the variance of the input. This suggests that the joint response of a cell population most 
accurately tracks input noise intensity. 



2 Model Description 

input 

f(V) + V2D$(t), ve(-/). (1) 



A single IF neuron with stochastic input is described by the Langevin equation: 

dV 
dt 



Here / defines the deterministic (drift) behavior, and <i;(f) a Gaussian stochastic processes with (£,(t)) = and 
(^(t)^(t')) = 8(t — t'). When the voltage reaches a threshold, V T , a spike is fired, and V is instantaneously reset 
to V R < V T . A spike may be followed by an absolute refractory period T, during which a neuron is insensitive to 
inputs, and V is held fixed at V R . 

This model can also be understood as a the diffusive limit of a population of cells receiving independent in- 
puts (Omurtag 2000). To model a pair of cells receiving correlated inputs, we assume that their membrane voltages V 
and W obey a pair of Langevin equations: 



V=f(V,W)+I v (t); I v (t)=!iv + V2D(VT^$v{t) + VcUt)) 



W=g(W,V)+I w {t); I w (t) = n w + V2D(VT^% w (t) + Vt%c(t)) 

The inputs, //(f), received by the cells are comprised of statistically independent stochastic processes E,v(t) and <^w{t), 
and a common input ^ c (t). The <!;,■ are again assumed to be Gaussian with (4,-(f)) = and (£<(*)£/(*')) = 8(t — t')8jj. 
The constant c, is the Pearson correlation coefficient between the inputs and lies between and 1 . For instance, for 
two leaky integrate-and-fire (LIF) neurons with common input, but no direct coupling, f(V, W) = — gv{V — V lest ) and 
g(V, W) = —gw(W — Wrest)- Each cell spikes when the voltage V crosses the threshold, V T and W T respectively. After 
each spike the voltage is reset to V R < V T {W R < W T for cell 2) , and is pinned to this value for the duration of the 
refractory period, Ty (Tiy for the second cell, see Fig . [TJ . For simplicity we will refer to the two as neuron V and W, 



although this can be understood as "populations V and W" (Harrison et al 2005 1. The joint probability density of the 
two voltages evolves on the domain D. = (V~°° \V T ) x (W~°° ,W' ). In theoretical studies it is frequently assumed that 
y-°° = = — oo. However, since we will be interested in numerical simulations, we assume that these quantities 
are large and negative. 

With U = (V,W) and F = (f,g) the Fokker-Planck equation corresponding to Eq. |2]) takes the form 

d t P(t,U)+6iv(F(U)P(t,U)-DMVP(t,U))=0, M = (3) 

forV e (y _0O ,y r )\y fl , W € (W~°°,W r )\V^. Here D is the diffusion coefficient and Mis the correlation matrix. This 
equation is coupled with reflecting boundary conditions at V = V ~°° or W = W~°°, and absorbing boundary conditions 
at both thresholds: 

(f(U) - Dd v - cDdw)P(t,U) | v=v ^ = 0, (g(U) -Ddw- cDdv)P(t,U)\ w=w _„ =0, 
P(t,U)\ v=V T=P(t,U)\ w=wT =0. 
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Figure 1: (Left) Domain of simulation. (Right) Circulation of probability mass through populations P, Ry, Rw and R. 



The presence of the refractory behavior in the IF model introduces additional complexity. If either neuron enters 
its refractory state, the corresponding voltage is fixed at the reset value, and the entire system effectively evolves 
according to a one-dimensional Fokker-Planck equation. During this time, it is possible that the second cell also 
crosses the threshold, fires and enters the refractory state. In this case the voltages are fixed at (V R ,W R ) and both 
neurons are insensitive to inputs until one of them exits the refractory state. 

To capture the behavior of neurons in the refractory period, we model the evolution of densities using three sepa- 
rate, communicating sub-populations, in addition to P(t,U) ( |Sirovich||2008[[Ly and Tranchina||2009| ): 

• Rv(t,r,V), the probability density of the fraction of the population in which only neuron W is in the refractory 
state, 

• Rw(t,s,W) the corresponding density in which only neuron V is in the refractory state, and 

• R(t,s,r), the density corresponding to both neurons in the refractory state. 

For all densities, t refers to the time since the beginning of the simulation, while r and s refer to relative times 
measured from the beginning of the refractory period for neuron V and W, respectively. Therefore, Ry (t,r, Vo)AVAr is 
the fraction of the population for which neuron W has been in the refractory period between r and r + Ar units of time, 
and the voltage of neuron V is between Vb and Vb + AV. The quantity R(t,s, r)AsAr is the fraction of the population in 



which neurons V and W have been in refractory periods for times in the intervals [5,5 + As] and [r, r + Ar 



The use of variables s and r is closely related to age-structured population dynamics models ( Iannelli 



1985). Indeed, s and r denote the "ages" of the refractory states for neurons 1 and 2 respectively, 
the circulation of probability mass between the different populations involved. 

Since the entire population is described by these densities we have for any time t 

V T f-W T rr w rW T 

/ P(t,V,W)dWdV+ / / R w {t,s,W)dWds 

v-°° Jw-°° Jo Jw-°° 

Ty rV T rt v rl w 



tneTEIJ l 



respectively. 



1994 Webb 



summarizes 



Jv 



rw r T w 

R v (t,s,V)dVds+ / R(t,s,r)dsdr= 1. (5) 
Jo Jo 



We next describe the evolution of the main population and the three refractory populations and how they are cou- 
pled to each other through boundary terms. The refractory populations evolve according to one-dimensional Fokker- 
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Planck equations, 



{d t + d r )R v {t,r,V) + div(f(U)R v (t,r,V)-DdvR v (t,r,V)) =0,re(0,%), 
(5, + d s )R w (t,s,W)+div(g(U)R w (t,s,W)-Dd w R w (t,s,W)) =0,se(0,z v ), 

(d t + d s + d r )R(t,s,r) = 0. 



(6) 



(7) 



In addition, mass from the main population P(t,U) is injected into populations Ry(t , r,V) and Rw(t,s,W) at r = 
and s — as neuron V and W cross threshold respectively. These source terms are described by 



R v (t,0,V) = -DdwP(t,V,W T ), R w (t,0,W) = -DdvP{t,V T ,W) . 



(8) 



While either neuron is in the refractory state, the other neuron can enter its own refractory state as well, providing 
boundary conditions to equation (j7]) for inward characteristics: 



R(t,s,0) = -Dd w R w {t,s,W T ), R(t,Q,r) = -Dd v R v (t,r,V T ) 



(9) 



Next, while both neurons are in the refractory state, neuron V or W may exit the refractory state while the other neuron 
remains in the refractory state. Using ] | v _ z : = lim r ^ z + | (z) — lim z ^ z - | (z) to denote the jump across point Zgl, 
we can express the contribution of population R(t,s,r) to the following source terms ( Melni kov] 1 1 993 [ |"Lindner| [200 1 ) 



[DdvRv(t,s,.j\ 



V=V R 



= R{t,z v ,s), [DdwR w (t,r,.)] 



W=W R 



■■R{t,r,x w ). 



(10) 



As either neuron exits the refractory period, it re-enters the main population modeled by the density P(t,U). This is 
captured by adding the source terms: 



\DdvP(t,.,V) 



V=V R 



:R w (t,T V ,W), [DdwP(t,V,.)] 



W=W R 



--Rv(t,Tw,V) 



The densities are also continuous across the reset potentials, so that 



[P(t,.,W)] 



v=v R 



[P(t,W,.j\ 



W=W R 



[R v (t,s,.)] 



V=V R 



[R w (t,s,.j\ 



0. 



w=w R 



Finally, reflecting and absorbing boundary conditions are imposed on Eq. ([6]), and (|7]i by requiring: 

(f(U)-Ddv)R v {t,s,.) lv=v -~ = (g(U)-Ddw)R w (t,s,.) lw=w -~=0 , 

R v {t,s,V T ) =R w (t,s,W T ) =0 . 
Notice that when there is no refractory period (Ty = Tw = 0), (|8]l-([TT|i reduce to the single boundary condition 

[DdvP(t,.,W)], v=vR = -Dd v P(t,V T ,V) , [Dd w P(t,V,.)], w=wR = -Dd w P(t,V,W T ) . 



(11) 

(12) 

(13) 
(14) 



3 Description of the Numerical Methods 

The numerical methods used in simulating the solutions of the model described in the previous section are not com- 
pletely standard. The anisotropy of the diffusion operator coupled with the absorbing boundary condition presents 
numerical challenges different from those encountered, for instance, when modeling phase oscillators (GalaiTeFal] 
2007). We therefore give a brief description of our approach here. 
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3.1 Finite Volume Method 



Three requirements in the numerical discretization of Eqs. (|3j)-(p"4|) are obtaining numerical probability densities which 
are accurate, nonnegative and integrate to 1. We dealt with similar difficulties in (Marpeau et al. 2009| l, and use an 



extension of that approach here. To obtain numerical densities which integrate to 1, we use conservative numerical 
schemes which ensure that the mass lost by a mesh-element is transmitted exactly to its neighboring elements. This 
ensures preservation of mass of the initial densities. Secondly, the drift operator is well known to be an obstacle when 
trying to combine accuracy, non-negativity, and stability of the numerical densities. Therefore, we use an operator 
splitting method (described below) that enables us to discretize the drift and diffusion operators separately. The 
discretization of the drift term is carried out by an upwind scheme whose accuracy is improved using flux limiters. 
However, the two-dimensional nature of the problem induces further difficulties: 

1 . Extra numerical diffusion is generated in the cross directions from the drift operator, leading to a loss of accuracy. 

2. Due to the correlation coefficient c, the diffusion matrix DM is anisotropic. The discretization of the cross 
derivatives dy W P commonly involves the inversion of matrices that are not unconditionally strongly diagonally 
dominant, which makes it difficult to obtain nonnegative numerical densities. 

3. Due to the refractory periods, the multiple boundary conditions (|6))-([T4|) drastically increase the algorithmic 
complexity, compared with the one-dimensional model in (Marpea u et al.j [2009| >. In particular, the densities Ry 
and Rw have to be obtained by discretizing Eq. (|6]l at each time step and age step. Moreover, the code has to 
gather all the phenomena and assemble the circulation between the four populations in an efficient manner. 

The technical of our approach are relegated to Appendix [A] and what follows is an outline. First, the time interval R + 
on which the solution will be approximated is partitioned into sub-in tervals (t n ,t n +\). We w ill denote the numerically 



obtained approximation of the population P at time t n by P". As in ( Marpeau et al.||2009| , P" +1 is obtained from P" 
by splitting equations into 

d t P + div(fP) =0, P\ v=v _ =P\ W=W _„=0, (15) 
d,P - div (DMVP) =0, 

(Dd v +CDd vw )P\ v=v _„ = (pd w +CDd vw )P\ w=w _» =P\ V=V - =P\ W=W ^=0. ' " 

The numerical solution is updated at time t n+ \ using 

p n+1 = y 2 (^i (/*»)), (17) 

where S?\ and S^i are approximation schemes for Eq. (jT3J and ([16} respectively, along with split interior conditions 
specified later. This technique allows us to develop specific numerical schemes which are adapted to each differential 
operator in Eq. Q. 

The numerical scheme is nonlinear explicit. A compromise between accuracy and stability is obtained by 
adding flux limiters to the upwind scheme. As discussed in the Appendix, we use the numerical scheme introduced 
in (Marpeau et al. 2009) in each direction, V and W. The time step is restricted by a Courant-Friedrichs-Lewy (CFL) 



condition (Courant 1928, Godlew ski and Raviart| [T990| , which provides stability and positivity preservation of the 



scheme by ensuring that the drift term does not shift the numerical solution by more than one mesh element per time 



step ( Godlewski and Raviart 1990 1. 



The scheme 5^2 is linear implicit. Using centered approximation of the derivatives dyyP, d^wP and a semi-center 
discretization of the cross derivative (see |A.2| >, the numerical solution is obtained by inverting a matrix that is strongly 
diagonally dominant as long as the mesh-size does not change too sharply between two elements. The scheme will 
remain stable and positivity preserving with any time step. The inversion of the linear system is carried out by an LU 
pre-conditioned gradient procedure. 

The one dimensional Fokker-Planck equations |6]) are discretized by using the one-dimensional scheme from 



(Marp eau et al.||2009| l (See A.4i. The main difference here is the presence of the age variables, s and r. Since the 
age evolves simultaneously with time, we just solve Eq. |6]i at each time step, regardless of age, and shift the age 
variable by one time step. Finally, Eq. ([7]) is solved exactly. The other boundary conditions given in Eqs. (|8]l-([T4]i are 
discretized as in one space dimension. The convergence criterion is satisfied as the residual of the numerical scheme 
decreases to a pre-defined value, 10~ 6 in our study. 
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3.2 Computing spike train statistics 

The statistics of the number of threshold crossings of an IF model are of special interest. Using neuroscience termi- 
nology, we will refer to each threshold crossing as a spike and the sequence of threshold crossings as a spike train. 
Let the stochastic set functions Ny(s,t) and Nw(s,t) denote the number spikes during the time interval [s,t] in cell 
V and W, respectively. The instantaneous firing rate of cell X = V, W is defined as the instantaneous rate at which 
the corresponding population of cells spikes at time t. In terms of the spiking probability of a single cell, this can be 
written as 

v x (t)= lim lpr(jVx(r,r + Ar)>0). 

Ar^O Af 

The conditional firing rate, v v \ w (x,t), is defined as the firing rate of cell V at time t + T given that W has spiked at 
time t, 

*Vlw(V) = lim -^-Pr (N v (t + x + At, t + x) >0\N w (t,t+At) > 0) 

Af^O Af 

and similarly for v W |y(T,f). The conditional firing rate can be normalized by the rates to obtain the spike train cross- 
covariance function 

C VW {X,t) = V W {t) (v V]W {x,t) - Vy (f + <)) , (18) 

which is a common measure of correlation between the activity of two neurons over time. In the study of neural 
coding, it is often useful to know the propensity of one cell to spike during some time interval after another cell has 
spiked. For this purpose, we define the conditional mean rate, 

1 f b 

S v \ w (a,b,t) = — r v v \ w (t,x)dx (19) 

(b-a) J a 

When the distribution of membrane potentials is in steady state, the spike trains are stationary and we can drop the 
explicit dependence on t to write Vx, Vyny(x), Cyw( T )> and S v i w (a,b) without ambiguity. 

Since action potentials are not explicitly modeled in the Fokker-Planck formalism described above, the spiking 
statistics must be calculated using properties of the probability density near threshold. The instantaneous firing rate of 
cell V can be obtained from the solution of the Fokker-Planck equation by taking the marginal flux over threshold, 

Vy{t) = -DP v {t,V)\ v=vT (20) 



where Pv(t,V) = J^P(t,V,W)dW is the marginal density of V. Thus, up to terms of 0((Al) ), the quantity Vy (f)Af 
is the probability mass that crossed threshold during the interval Af. Equivalently, it equals the probability that a cell 
fires during this interval. The instantaneous firing rate of W is defined analogously. 

The conditional firing rates are obtained by first calculating the conditional flux immediately after a spike in cell 
W at time t , 

■W(f,V) := -Dd w P(t,V,W)\ w=wT . (21) 

This conditional flux is then normalized to give the conditional density, J CO nd(t,V) J C ond(t,V) j fy-*, J con a(t ,x)dx 
and used as an initial condition for the 1 -dimensional Fokker-Planck equation, 

d T P, (V, x) = —dv(f(V) - D dvP { (V, x) ) . (22) 

As the solution of this equation evolves, the conditional firing rate of V is given by Vy\w(x,t) = —DdyP\ (V, x). The 
conditional firing rate can then be normalized, cf. Eq. ( [IS) , to get the cross-covariance function, or integrated, cf. 
Eq. ( [T9| ), to get the conditional expected mean rate. We experienced convergence problems with the derivative of the 
finite volume solutions at the upper corner, (V T ,W T ), of the spatial domain. Due to the convergence issues discussed 
in Appendix|A.6| the finite volume approximation to the conditional firing rate does not converge when x is very small. 
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Figure 2: Stationary probability densities for a pair of LIF neurons: Results from finite volume simulations (first), 
MC simulations (second), the difference between the two approximations (third column). To test convergence of the 
finite volume method, we used a coarse (100 x 100 elements in the unit square), and a fine grid (200 x 200 elements). 
The fourth column shows the L 1 norm of the difference between the equilibrium distributions obtained using the finite 
volume and MC simulations. Parameters: for the top panel, jXw = Hv = 0.5, D = 0.05, c = 0.5, T = 0.5; for the middle 
panel, p w = /iy = 0.5, D = 0.05, c = 0.9, T = 0.5; for the bottom panel, pL V = 1.2, /iw = 0.6, D = 0.05, c = 0.3, 
T = 0.2. The first three columns were obtained using a coarse (100 x 100) grid. 



4 Validation of the numerical solution 



As the finite volume numerical scheme we developed is novel, we first compare its output to that obtained using Monte 



Carlo (MC) simulations (See Appendix A. 5 I. We consider both stationary and non-stationary inputs. 

As an example we choose the case of two LI F neurons which corresponds to setting f(V,W) = — V + fly, and 
g(V,W) = -W + Hw in Eq. ^ ( Burkitt 2006 1. When fly < V T and )i w < W T the cells are in the fluctuation 
dominated regime, and firing is due to large excursions of membrane voltages from the mean. As shown in top and 
middle panels of Fig. [2] the finite volume method provides an excellent approximation of the stationary distribution 
when the input to the two cells is constant in time. As the correlation between the inputs to the two cells, measured 
by c, increases, the membrane potentials become more correlated, and their joint probability density is stretched along 
the diagonal. 

When fly > V T or jiw > W T it is the DC component of the input current that drives the cells over threshold. This 
situation is somewhat more challenging to simulate, since much of the mass of the invariant distribution lies close to 
the threshold. The gradient of the solution close to the boundary becomes large. Together with the Dirichlet boundary 
conditions, this causes larger errors in the numerical approximation close to the boundary. The bottom panels of Fig. [2] 
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show that the finite volume method still performs well in this situation. 

We can change the drift term in Eq. Q to simulate a different integrate-and-fire model. In particular, the quadratic 
integrate and fire (QIF) model is obtained by setting f{V,W) = V 2 + Hv,g(V, W) = W 2 + jXw ( fErmentrout and Kopell| 
1986 Brunei and Latham 2003 1. Fig. [3] demonstrates that the finite volume numerical scheme performs well in 
computing the invariant distribution for this model. 

The finite volume scheme was designed to compute the evolution of the joint probability density of the two sub- 
threshold voltages in time. Stationary distributions were presented here for ease of visualization. A comparison of time 
dependent solutions obtained using finite volume and MC methods is available online at |http : / / www . math . uh . | 
edu/~ j os ic/myweb/ research/papers /FV/| The animation shows the time dependent density from t = 100 
to t = 1 10 for a pair of LIFs with pi = |sin(f)|, c= |sin(f)|/2, andD = 0.1|sin(f)|. 



Finite Volume 



Monte Carlo 




Figure 3: Stationary probability densities for a pair of QIF neurons: finite volume method (left), MC simulations 
(right) when /i v = Hw — —0.1, D = 0.1 and c = 0.3. 



5 Gaussian approximation 



The LIF model is ubiquitous in stochastic modeling of excitable systems primarily due to its mathematical tractability. 
Closed form expressions have been obtained for the stationary density, spiking statistics, and linear response properties 



of the one neuron model (Lindner 2001 1. However, exact closed form expressions are not known for the two neuron 
model with c ^ discussed in Section[2] The numerical methods we describe here offer a way of exploring the behavior 
of the LIF model in the absence of analytic solutions. However, even with fast numerical methods, exploring large 
regions of parameter space may not be possible. Approximate analytic solution, are therefore frequently necessary to 
gain a deeper understanding of the model. 

Much recent work has focused on deriving such approximations using perturbative methods. Linear response 



theory was used to study the dependencies in the output of a call pair receiving correlated input ( Lindner 2001 Ostojic 



|et al.||20"0"9"t|de la Rocha et aT| 2007 ; Sh ea-Brown et al.||2008| l. These solutions involve integrals that must be evaluated 
numerically. Simpler approximations can be obtained by ignoring the threshold and reset boundary conditions when 
the neurons are in the fluctuation dominated regime, and firing rates are low. Neurons in the cortex may reside in 
this regime under typical conditions (Ringach and Malone 2QQ7\. Previous approximations obtained in this regime 



required smoothness assumptions on the trajectories of the membrane potentials (Burak et al. 2009 Tchumatchenko| 
et al. 2Q1Q| >. Since solutions to Eq. Q are nowhere differentiable when D > 0, a different approach must be used for 
the LIF driven by white noise inputs. We next describe such an approximation. (We note that a similar approach has 
been used to examine the response of integrate-and-fire neurons driven by filtered Gaussian noise (Ba del et al.||2010) .) 

When firing rates are small, the boundary conditions have a small impact on the solution of Eq. Q and an ap- 
proximate solution can be obtained by solving the free boundary problem (V T ,W T —> °°). Since firing is rare, the 
amount of time spent in the refractory states is negligible and the refractory period can be ign ored (T = 0). Un der this 
approximation, the stochastic process (V(t),W(t)) is an Ornstein-Uhlenbeck process in K 2 (Gardiner 1985 I. Given 
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bivariate Gaussian initial conditions, the solution to the Fokker-Planck equation at any time is a bivariate Gaussian and 
can be computed in closed form. This Gaussian approximation is accurate when X T — jJ-x^> \/2Z) for X — V,W. 

For simplicity, in this section we assume that the two neurons receive statistically identical inputs so that fly = 
j±w = M- We further assume that the neurons are dynamically identical so that gy = gw, V mst — W res t, and V T = W T . 
The analysis is similar in the asymmetric case. Without loss of generality, we rescale space so that V res t = W lest = and 
V T = W T = 1. To simplify calculations, we also time in units of the membrane time constants so that gy — gw = 1. 

For instance, the marginal or conditional firing rates can be approximated by the flux of the time dependent Gaus- 
sian distribution over threshold. As shown in Appendix [B] this flux can be written in terms of the mean (m) and the 
variance (a 2 ) of the Gaussian, and the input diffusion coefficient (D) as 

W.D):-^.^. (23, 



The steady state firing rate, Voo is obtained by taking m = and a 2 = D to get Voo = J(jJ,,D,D) = -j=e a where 



2007 



Lindner 



a := (1 — /i) / %/2D- This expression can also be obtained using large deviation methods (Kampen 
[200T]|Shea-Brown et al.|[20Q8] >. 

From an approximation of the conditional firing rate, the cross-covariance function can be obtained as (see Ap 
pendix|B]), 

1 2 „ 2 / e'-^T^ 



C vlw (T) = v^H(T)-v^ = -a z e 



7t I Vl -c 2 e- 2r (c + e T ) 

In Fig. |4] we compare this approximation to the cross-covariance function to the cross-covariance function obtained 
from finite volume simulations. As expected, we find that the two agree well when firing rates and correlations are 
small, but disagree when jx, D or c are larger. 

Further expressions for other stationary and non-stationary spiking statistics under the Gaussian approximation are 
derived in Appendix[B] We use these approximations to examine the response of a pair of cells to time-varying inputs 
next. 



5.1 Response to step changes in the input - single cell response 

If a pair of cells responds rapidly to a change in an input parameter, then the output of the cell pair can accurately 
capture the information present in a time-varying input signal ( Silberber g et al.j [2004 ; Masuda, 2006). It is therefore 
useful to understand how the spiking statistics of a neuron or pair of neurons respond to changes in input parameters, 
jJ.,D, and c in Eq. 0. The response of the cell pair is measured by their joint firing rate (vy, Vw), and we first examine 
how rapidly this response can track changes in the inputs to the model. 

We start by examining the response of a single cell. Fig. [5] Left shows the time dependent firing rate after a 
step change in the mean, (ji, light line), and variance, (D, heavy line). We compare the response for the Gaussian 
approximation, derived in closed form in Appendix[Bj to the result from finite volume simulations. After a change in 
parameters, the distribution of (V,W), and therefore the spiking statistics, relax exponentially to a new steady state. 
However, the speed of this relaxation depends on the parameters that are changed. After a step change in the mean 
input, ji, the firing rate relaxes to a new steady state with a time constant of 1 (i.e., one membrane time constant). 
After a change in the variance, D, the firing rate jumps discontinuously, then approaches the new steady state value 
with a faster time constant of 1/2 (see Appendix |B|i. The fact that changes in the variance of the input are tracked 



faster than changes in input intensity is a fundamental property of the Ornstein-Uhlenbeck process (Gardiner 1985). 
Therefore, changes in variance can be tracked more faithfully than changes in the DC component of the input. Related 
observations are made in ( Khorsand and Chance , 2008 1 Hasegawa , 2009 ; Silberberg et al. , 2004 ; Masuda , 2006 ), where 



mainly the discontinuous change in output firing rate in response to a step change in input variance was examined. 

It also follows that a transient pulse change in D results in a larger transient in the firing rate than a comparable 
pulse change in ji. This prediction is verified for the Gaussian approximation and for finite volume simulations in 
Fig.|5|*tofa. 
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Relative L 1 difference 




0.25 



Figure 4: Left: Cross-covariance functions for pi = 0.1 and fi = 0.25 when D = 0.05 and c = 0.2. The solid lines were 
obtained from the Gaussian approximation and the dashed lines from finite volume simulations. Right: The relative L 1 
difference between the Gaussian and finite volume cross-covariance functions (L 1 difference divided by the L 1 norm 
of the finite volume result) for jti e [0.075,0.25] and c € [0.075,0.3 1. Th e L 1 norm was computed for T > 0.15, due to 
the convergence issues of V v \w for small T discussed in Appendix A. 6 The cross-covariance function has units Hz 2 



In this figure and in Figs. [5] and [6] the axes are labeled assuming a membrane time constant of l/gv = 1 /gw = 5ms 



v(t) (Hz) y( t ) (Hz) 




t (ms) 



Figure 5: Left: Instantaneous firing rate after a step change in input statistics. Result from the Gaussian approximation 
are plotted as solid lines and finite volume simulations as dashed lines. For time t < 0, the parameters were set to 
ji = and D = 0.03. At time t = the parameters were changed to = 0. 1334 (grey) or D = 0.04 (black). The values 
of ji and D were chosen so that the steady state firing rate after the change in parameters was the same whether or D 
was changed. Right: Instantaneous firing rate after a pulse change in input statistics. Same as Left, but the parameters 
were changed back to fi = and D = 0.03 at time t = 2.5ms. Note that a step change in the input variance, D, 
results in an instantaneous jump in the firing rate, followed by a continuous relaxation to the steady state. To illustrate 
the quantitative accuracy of the Gaussian approximation, we used parameters that resulted in low firing rates. Fig. [6] 
illustrates that, while the Gaussian approximation is less accurate when firing rates are moderate, the approximation 
can still capture the qualitative behavior of the spiking statistics. 
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5.2 Response to step changes in the input - joint response 



We next examine how joint response of the cell pair in response to a step change in the input. The cross-covariance 
function, defined in Eq. ( fT8) l, is commonly used to measure dependencies between two spike trains over time. Fig. 
[4] compares the Gaussian approximation of the cross-covariance function to that obtained using finite volume simu- 
lations. As expected, the two results agree when firing rates are low. However, as fi or D increase, the firing rate 
increases and the Gaussian approximation breaks down. 

Fig. [6] Top shows the two-point conditional firing rate after a step change in the parameter D. This function 
completely characterizes the second order correlations of the two cells over time. Such a plot would be computationally 
prohibitively expensive to obtain using direct Monte-Carlo simulations, especially when firing rates are low. The 
parameters for Fig. [6] were chosen so that the Gaussian approximation does not agree quantitatively with the finite 
volume simulations. However, as shown in Fig. [6] Bottom, the Gaussian approximation successfully predicts the 
qualitative behavior of the bivariate spiking statistics with changing input parameters. 

A pulse change in D has a larger impact on the propensity of the cells to fire together than a comparable pulse 
change in jj, (compare to Fig. |5]/?/g/2f). A pulse change in c has an intermediate impact. If a downstream cell that 
receives inputs from cells V and W is sensitive to synchrony in its inputs (Salinas and Sejnowski 2000), then the cell 
would response more quickly and strongly to changes in D or c than to changes in jj,. 

This suggests that downstream cells that act as coincidence detectors are most sensitive to upstream changes in 
input variance, and respond more weakly to changes in input intensity. 



6 Discussion 



Population density methods have a long history in neuroscience. They have been used to study both the statistics of 
responses of single neuron (Tuckwell 1988 ), and neural populations (Harrison et al 2005) . Studying the evolution of 
the ensembles, rather than tracking individual neurons has several advantages: While the dynamics of each individual 
cell is stochastic, their probability density evolves deterministically. The probability density is therefore easier to study 



analytically and numerically. For instance, both linear response methods ( |Lindner 2001 de la Rocha et al. 



Ostojic et al.| |2009[>, and the Gaussian approximation discussed here (see also (Burak et al. 2009; Tchumatchenko 
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et al. 2010 |Badel et~aL 2010) )), are obtained by considering the evolution of large populations in the diffusive limit. 
Simulating the evolution of population densities is typically orders of magnitudes faster than simulating the evolution 
of each individual cell in a population (Nykamp and Tranchina 2000 1 Omurtag 2000 ). For instance, obtaining the 
two-point time dependent firing rate, v CO nd(T,f) shown in Fig.|6| would not be feasible using Monte Carlo methods on 
an average machine today. 

We have concentrated on a simple version of the model to keep the presentation relatively concise. For instance, 
it is easy to consider cells receiving different, potentially time dependent drives, fly (t) , fiw (t) , Dy (t) and Dw(t). We 
have mainly considered drift terms of the form, f(V,W) = f(V), and g(W,V) — g(W), in Eq. Q, so that the two 
populations were uncoupled. Coupled populations have been considered earlier (Nyk amp and Tranchina||2~000| l. We 
hence concentrated on examining the effects of anisotropic diffusion. However, the numerical methods we described 
can easily handle coupling between cells in the sub-threshold regime. This could be used to examine the interplay of 
correlated inputs and cell coupling ( Schneider et al. 2006 ). However, we do not know whether there is a direct way to 



include super-threshold coupling in the present diffusive approximation. We also note that high firing rates can lead to 
steep gradients of the probability density close to the boundary and convergence problems in the numerical methods 
we used. This suggests that numerical techniques will have to be developed further to accurately capture the response 
of IF neurons in this regime. 

We have used our numerical and analytical approach to examine the best way to transmit information in a pair 
of cells. We found that both the single cell response and the joint response tracks changes in noise intensity more 
accurately than changes in the mean drive or correlations in the cell inputs. Thus input variance appears to provide the 
best channel to code information at the single cell and population level. 

The numerical methods we have developed can be used to further examine how the output of a cell pair reflects 
their interactions and dependencies in their inputs. It is of particular interest how cells respond to signals that vary in 
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Figure 6: Top row: The two-point time dependent firing rate, V con( j(T,f ) after a step change in D with D = 0.1 for t < 
and D — 0.2 for t > 0. Parameters ji = and c — 0.1 were held constant. Values for T < 0.5ms were omitted due to 



the convergence issues discussed in Appendix A.6 Bottom row: The conditional mean rate, S, after a pulse change in 
parameters. The parameters changed from D,c,fi — 0.1,0.1,0 for t < to D = 0.2 (black), c = 0.435 (dark grey), or 
ji = 0.293 (light grey) for t 6 [0,2.5]ms, then back to D,c,ji = 0.1,0.1,0 for t > 2.5ms. The conditional mean rate, 
calculated according to Eq. fL9] l with a = 0.15ms and b = 0.5ms, measures the propensity for cell V to spike within 
the first 0.5ms after cell W spikes. We chose a = 0.15ms to circumvent the convergence issues discussed in Appendix 



A.6 Left column: finite volume simulations. Right column: Gaussian approximation. 
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Figure 7: A schematic depiction of the subdivision of the domain (V °°,V T ) into subintervals , V i+ i ). The reset 
voltage V s is at V, w . Similar notations are used for the second neuron voltage W. 

time. The finite volume method we described is well suited to this task, as it is designed to capture the time-dependent 
response of a cell pair. 

A Numerical schemes 

In the following we provide a description of the finite volume method used in the numerical simulations. Some of the 
details of the implementation are not standard. As we are not aware of a similar treatment of this type of equation, we 
give a detailed discussion of the novel aspects of the algorithm. 

The time steps are defined by At„ = t„ + \ — 1„. When there is no ambiguity, the time step is denoted by Af. The 
intervals (V~°°, V r ) and (W~°°,W T ) are partitioned into Ny and Nw sub-intervals respectively. We denote AV; as the 
i th step in the V -direction and AWj as the /* step in the W-direction, for i = 1, . . . ,Afy, 7 = 1,... ,Nw- Our quadrilateral 
mesh elements, Q, -j, are then defined by 

Qu = (v i ^,v iH )x(w jA ,w j+k ), 

with V._ i = J AVi, Wj_ i = Y^jZl AW/. Our mesh-points (Vi,Wj) are the centroids of the cells, thus V; = V f _ i + \ AVi, 
Wj = Wj_i + \Wj. We make sure that there exist two indices ig and jx such that Vj R = V R and Wi R = W R , which 

means that the two reset potentials fall exactly on some mesh-points, see Fig. [7] For every function £ defined on 
(0,r) x Q., the notation £j£ ^ stands for the approximation of t, (t n , (V a , Wjg)), for Of = i, z'± J3 = We denote 

t," as the sequence {§",•};,;. 

A.l Treatment of the drift operator: scheme S?\ 



The advection equation in ( 15 1 is discretized by using the one-dimensional numerical fluxes in (Marpeau et al. 2009) 
in each direction. We set 

p "t 1 =^--t L (<' r<" 1 P-^HV (24) 

'» ' J Ay; 1 2J Awy ! J+2 'J J 

where the numerical fluxes are defined by 

2 Av ;+ i V ! +zJ ,+ i>} AVi l +i>J l+ 2'J Av,-+i / 
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for 1 < i < N v - 1, 1 < j < N w and 



, — o + P" . -I- CT~ p." . , 



1 ^"i ';+! / ^^V 1 Aw., i \ 

9 (<;+ ' Aw M r n + ' > V^) - *n + ' A ^ +1 ^ ?/+ ' ' T^)) ' (26) 

2Aw ( j + i V ij+2 v Awj ' >>J+2 ,J+ 2 Aw,+i ' J 



for 1 < z < M,, 1 < / < AW - 1, with notations AP" , := P" , - if,, AP" , := P/ 1 , , , - /f ,, 

f + , AP. i . / " . AP.^ i . g + | AP. . i g~ , AP. . i 

p = J i-j.j >-l'J m = J i+ l 1 .j >+l>l p _ 6 iJ-j •'J-l m _ 6 iJ+j 'J+2 

r i+iJ~f + , AP.^i .' i-kJ r. AP. , . ' r 'J+z ~ g + | AP. ., i ' r, ''./-2 _ p~ ,AP . i ' 

where the limiter function <p is defined by (p(a,b) = 2b max (0,min(l,2a),min(a,2)). 

To comply with Dirichlet boundary conditions in (JT5|, we further impose stftj = + , . = J2^" = + , = 0, 

for 1 < i<N v , 1 <j<N w . 



Proposition 1 The numerical scheme ( 24 1 is non-negativity preserving under the CFL condition 

Av,- Awj \ Avi J V Awj 

The proof of this proposition is similar to the proof of the corresponding one-dimensional result in ( jMarpeau et| 
|aL)|2009| ), and we therefore omit it here. 

A.2 Treatment of the diffusion operator: scheme S?i 

Our approximation of the solutions to ( [T6| in two space dimensions is given by the implicit scheme 

Av^y „+l _ ( ^+, _ ) _ _ ) = A ^ P » + 5i iR Sf" +l + 8j jR Sj' n+1 , (27) 

At '>■> >+2-] i— j j v «,;+2 'J-z Af J '* J • / '' / * 1 

where 8k { .k 2 is the symbol of Kronecker and sj'" +l and sj'" +l account for the re-injection condition |TT[ i, see subsec- 



tion 



A. 3 The numerical diffusive fluxes are defined by 



pn+l _ pn+l pn+l _ pn+l pn+l _ pn+l 

=£) fmj ' ) . (28 ) 

<+ 2 J Av,.,i 2 V Aw., i Aw. i 

'+2 J+J J—J 



for 1 < i < N v - 1, 1 < j < N w and 



pn+l _ p«+l _ p«+l _ pn+l pn+l _ pn+l 

an+1 =D lU+l lid I C ( J+hi+1 •■J+ l I lthi "\ (29) 

Aw,, i 2 V Av., i Av, i ' ' 

J+j i+ 2 I- j 



for 1 < i < N v , 1 < j < N w - 1. 



Remark 1 Notice that the second term of the right-hand side in (28 I stand for a centered finite difference discretiza- 
tion of the cross derivative CDd^ w P on the right vertical interface of Qjj. Other numerical schemes have been im- 
plemented in (Bruneau et al. 2005 ' Rasetarinera. 1995' Bourgeat and Kern. 2004), but yield unconditionally positive 
off-diagonal coefficients in the diffusion matrix, therefore producing negative undershoot near sharp solution gradi- 
ents. When the neurons are strongly correlated (i.e.: when C ~ \), the gradients of the solution can be very sharp. 
The advantage of our method is that all the off-diagonal coefficients are nonnegative where the mesh is uniform, which 
means, the region where the solution is not in practice. Using a similar remark in ( 29 ), the resulting numerical 
scheme ((27]) is nonnegative in realistic applications. 
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Figure 8: A schematic depiction of the time dependent subdivision of the domain (0, Ty) into subintervals 



The boundary conditions in (116} are implemented as 

7,'J 



23' 



•n+1 

1 



fl+\ 



NvJ+2 



21+1 



D- 



pn+1 



AV/2 



D- 



D 



pn+l 
i+\,N w 

AV. 



pn+l 



Tfi+1 



CDfO-P, 



,«+l 
i+lJVw 



= 0and 

A _ pn+l 

AW N +1 

N w +2 

pn+l 



pn+l 



2 \ AW Nw /2 



= D 



pn+l _ 
Ny,j+l 1 Ny,j 



□n+1 
'Ny.j+l 

2 V AV Nv /2 



pn+l 



pn+l 
r AV-lj 



A.3 Treatment of the re-injection condition ( [XT) 

Since the time variable t and the age variable s G (0, Ti) evolve together, the domain (0, Tj ) is dynamically partitioned 
into sub-intervals (^,^+i) sucn tnat ^ = an d -^+1 = rnin(sjj +ASj,Ti), where the age steps As", match the time 
steps as follows: for all n, As" = At„ and As" k = Af„_£, k = l,...,n— 1. We set = max{£, < Tj}, see figure [8] 
In the same way, we discretize (0, T2) into sub-intervals (rj|!,r£ +1 ) such that r" k =Q and rjj, j = min(r^ + Ar k , Ti), with 
Ar" = Af„, ArJ! = Af„_fc, = 1, . . . ,n - 1. We set = max{£, r£ < T2}. 
Then, defining the piecewise constant functions 

Ri,i(t) = + iK- K X (tn , n+i] {t) , R 2J (t) = f^f K % n , tn+l] (t) , 

n=0 n=0 

the quantities 

Sf := / Ri ti (t)dt and Sj :=[ * R 2 ,j(t)dt 

are injected into p7j ). 



A.4 Treatment of the one neuron boundary condition (|6]>— (|T0]> 

Equation ^ is discretized with the one-dimensional numerical scheme in (Marpeau et al. 2009 1. We use the operator 
splitting technique 



R 



2^+2 _ ftn.k _ 



At 



K+>-<.) 



and R 



n+l,i+l 
1,1 



= R 



+ S i;iR S 



W,n+1 



(30) 
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where the advective and diffusive numerical fluxes are 



R 



n.k 



1 



2 AV i + i 



(f + , AVj(p(r p 
\ J i+l,N w ,r V f+ 



:■/•/ 



.n+l 
t+i 



D n+1.*+1 D n+l,i+l 
K l,i+l ~ K U 



for z = l,Ny — 1, and , 



0, , 



AV+ 



Av 



'+2 
D _ftii+l,k+l 



(31) 



We use similar arguments to discretize Eq. (|6|). Then, since the age steps As? match the time steps A?„_£, equation 
(|7|i is solved exactly as — R"^- 1 , and its boundary conditions (9 1 translate into R n <°< 1 



R" 



k.O . 



D 



R 



n.k 



2AV> 



D R nJ 



N V 



1,AV' 



2AW Nw "2,N W 



Finally, the first re-injection condition in ( 10 1 is taken into account by setting S W ' n+l ' 1 := R n+1 > KS > 1 in (|30|>. The 
second re-injection condition in ( [T0| is discretized in a similar fashion. 



A.5 Comparison with solutions obtained using Monte Carlo methods 

Monte Carlo simulat ions were performed by inte grating the Langevin equations Q using the Euler-Maruyama method 
with time step 10~ 3 ( |Kloeden and Platen 1992 1. Histograms in the V — W plane were typically created using 2 x 10 7 
points sampled from along, simulated trajectory and using the same grid on the domain £2 = (V~°°,V T ) x (W~°°,W T ) 
as in the finite volume simulation. To compare the results of the finite volume and Monte Carlo methods, we computed 
the L 1 norm of their difference by taking the absolute value of the difference at each grid point and summing over the 
domain £2. 



A.6 Calculating spike train statistics from finite volume simulations 

To obtain spike train statistics from finite volume simulations, we used the definitions from Sec. |3.2| We solved the 
two-dimensional Fokker-Planck equation using the numerical scheme described above. To solve the one-dimensional 
Fokker-Planck equation (for example, to obtain instantaneous or conditional firing rates cf. Eq. <[22]>), we used the 
one-dimensional scheme described in (Marpeau et al. 2009[ > or we used the two-dimensional scheme and calculated 



the marginal density from the two-dimensional density as Pi (t, V) — f^-°° P(t ,V,w)dw. 

However, we met with convergence problems in calculating the conditional flux J CO nd(t , V) using Eq. ( f2"Tj ). These 
are due to the fact that when W is close to the threshold value W T , and the two neurons fire nearly synchronously, 
the re-injection process in the V and W directions is relatively complicated. Given these convergence issues, we 
ignored the first 0.5ms in the left panel of Fig. |4]and 0.75ms in right panel of Fig. |4]and Fig. |6]when we computed the 
conditional firing rate Vyi w (T,f). 



B Gaussian approximation of the LIF 

In this section we derive an approximation of the spiking statistics for the LIF in low firing rate regimes. Recall that 
the LIF is defined by taking f(V,W) = -g V (V - V mst ) and g(V,W) = -gw(W - W mst ) in Eq. g). For simplicity, in 
this section we assume that the two neurons receive statistically identical inputs so that fly = Mw = M- We further 
assume that the neurons are dynamically identical so that gy = gw, V re st = Wrests and V T = W T . The analysis is similar 
in the asymmetric case. Without loss of generality, we rescale space so that V res t = W KS t = and V T = W T = 1. To 
simplify calculations, we also time in units of the membrane time constants so that gy = gw = 1. 

When a := (1 — fx) /v2Dis large, firing rates are low and the boundary conditions at threshold have a small impact 
on the distribution. Also, since the cells only spike rarely, the refractory period can be ignored. In such regimes, the 
solution of the full problem is approximated by the solution of the free boundary problem. This approximation, which 
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we call the Gaussian approximation, is accurate in the limit a — » °°. In this case the membrane potentials (V(t),W(t)) 
are described by an Ornstein-Uhlenbeck process on R 2 . Such processes are well-understood and the spiking statistics 
can be computed exactly, as we show below. 

Assume that the initial distribution P(Q, U) is abivariate Gaussian with marginal means m(Q) =E[V(Q)} =E[W(Q)}, 
variance <7 2 (0) = var(V(0)) = var(W(0)), and covariance y(0) = cov(V (0), W(0)). Then the solution at any time t > 
(in the absence of boundary conditions) is Gaussian with mean, variance and covariance given respectively by 



m(t) 


= e-'m(0) + 


(1 — e f )m(°°), 


a 2 (t) 


= e- 2 'a 2 (0) 


+ (l-<r 2r )a 2 H, and 


7(0 


= e- 2 ' y(0)4 


(i- e - 2 ')rH 



where 

m(°°) = /I, (7 2 (°°) = D. and 7(00) =cD 
are the steady state mean, variance, and covariance. 



B.l Conditional firing rate and spike count correlation in the steady state 

The results above can be used to derive an approximation of the steady state conditional firing rates. Since the joint 
distribution of (V,W) is a bivariate Gaussian, the distribution of V given that W = W T = 1 is a univariate Gaussian. 
The conditional mean and variance are m c (Q) — c(l — ju) + /I and (7 2 (0) = D(l — c 2 ) respectively. As time evolves, 
the conditional density of V relaxes to its steady state. The density during this relaxation is a univariate Gaussian with 
mean 

m c (T) = e~ T m c (0) + (1 - e~ x ) m c (oo) 

and variance 

a 2 (T) = e - 2T a 2 (0) + (l-e- 2T KH 

where m c (°°) = [l and G 2 {°°) = D are the stationary mean and variance. Note that for c near 1, m c (0) is near V T = 1 
which violates the assumptions of the Gaussian approximation, namely that the mass near threshold is small. In this 
case, the conditional flux across threshold is large even if the marginal fluxes across threshold are small. Thus, for the 
approximation of the conditional firing rate to be accurate, we must assume that a is large and that c is small. 
The conditional firing rate is simply 

v v \w{t) =^(A-(tW(t),D) 

where J(jj,,a 2 7 D) is as defined in ( |23j ). This expression does not depend on t because we have assumed that the 
two-dimensional distribution is in its steady state. Later, we look at the two point conditional firing rate outside of the 
steady state. 

The steady state cross-covariance function is obtained from the conditional firing rate cf. Eq. ([18) to obtain 



1 / e 

CW( T ) = Voo(#(t) - Voo) = -a 2 e~ a ~ 



To first order in c this gives, 
The asymptotic spike count correlation, defined by 



Vl-c 2 e- 2z (c + e T ) 



R(x) = -a 2 (2a z -l)e- 2a - r + o(c 2 ). (32) 



= cov(N v (t),N w (t)) 
'->» y/var{N v (t))var(N w {t)) '' 
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can be written in terms of the conditional firing rate as (Shea-Br own et al!)|2008| ) 

_ 2/ °°(vyny(T)-V,)rfT 

P ~ CV 2 

where CV is the coefficient of variation of the spike train inter-spike intervals. In the low firing rate, a 
CV — > 1 and therefore, to first order in c and Voo, 



»2/ (v V | W (t)-v« 
Jo 

= -^=2a(2a 2 -\)e 



\dx 



limit, 



(33) 



Since both Voo and the correlation susceptibility, T := are functions of the single parameter a and since Voo is 
monotonic with a, we may conclude that T is a function of Voo to first order in c and v^. This same conclusion was 
reached in Shea-Brown et al. ( 2008| l using linear response theory, though the expression derived for p, 

2 



a 2a- 



, e 



(34) 



dp 



• 4ca as a — > °° (i.e., as Voo — > 0). Comparing these two approxima- 



differs from Eq. ( [33) . For both expressions, 

tions to a more accurate linear response approximation (also from (Shea-Brown et al. 2008 )), we found that Eq. p4| ) 
is more accurate than Eq. &3\. 



B.2 Time dependent input statistics 

We will now investigate how the spiking statistics track time dependent changes in the inputs. When the input param- 
eters to the neurons are time dependent, the two dimensional density p(t,V,W) at any time t is a bivariate Gaussian 
whenever the initial condition is a bivariate Gaussian. Thus we can use the same methods as above to derive the 
time dependent spiking statistics. To illustrate the effects of time-dependent inputs, we concentrate on a simple time- 
dependent input model. We assume that each cell receives input with mean /Iq, diffusion Do, and correlation cq for 
t < and, at time t = 0, the input parameters change instantaneously to D\, and c\. At some later time to > 0, 
the inputs change back to the original values, /Zo, Dq, and cq. A small value of to models a pulse change in the in- 
puts. Taking t q = °° models a step change. The discussion here can easily be generalized to arbitrary time-dependent 
input (e.g., sinusoidally varying inputs) by solving a simple linear ODE for the time dependent mean, variance, and 
covariance ( |Gardiner||1985| l. 

We assume that for time t < 0, the distribution is in its steady state so that 

m(t) 

t<0 




At time t = 0, the input statistics change and the mean and covariance matrix begin to track this change. In particular, 

m(t) = e-'ii + (l-e-')ii u "I 

d 2 (f) = e - 2 'D + (i- e - 2 ')£> 1 I te[0,t ] 

7(f) = e - 2 ' C0 A) + (l-e- 2 ')ci£>i J 

At time to, the input statistics change back to jj,o, Do, and cq and the distribution relaxes back to its steady state. In 
particular, 



m(t) = e-(' - *>)m(fo) + f 1 - e~ ( '~' o) ) Ho, 
a 2 (t) = e - 2 ( r -'o)a 2 (ro) + (l-e- 2 ('-'°>)z)o \ t>t 



7 ( f ) = e - 2 ('-'o) 7 ( fo ) + (i _ e -2(t-t )^ CqDq 
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where /i(?o), (j2 ( f o) and 7( f o) are given by the previous set of equations. The variance and covariance of the solutions 
change with a time constant that is twice as fast as the time constant with which mean changes. This is a known 
property of Ornstein-Uhlenbeck processes (Gardiner 1985). 



B.2.1 The time-dependent firing rate 



We now investigate how the firing rate changes in response to a pulse or a step change in the input statistics. The firing 
rate at time t is given by v(t) — J(p(t),0 2 (t),D t ) where J(ji, <7 2 ,D) is defined in ( |23] l, p,(t) and G 2 (t) are as derived 
above and 

(d t<0 
D, = \d x te[0,f ] 
[Do t > t Q 

is the time dependent diffusion coefficient. 
We can simplify the expression to get 



Vy(f) 



( "(0 „-a 2 (t) 



a{t)- a Ht\ 



e- 2, D +( 1 -«~ 2 ') D1 J ^ n 
Do 



K I (l+e- 2 '-e- 2 ('-'o))£> Q +( e - 2 ('-'o)- e -2r) Di y sfn 



«W.-a 2 (f) 



t <0 
te [0,t ] 

t>to 



where 



a{t) 



l-m(t) 



multiplying the -j+e 
a factor of 



a(t)„-aHt) 



Note that a(t) changes continuously with t. Thus any discontinuities in the expression above are from the factors 

term. In particular, the firing rate jumps discontinuously by a factor of g 1 at time and by 
at time to- If we change the mean of the input signal, but do not change the variance of 



the input signal (by setting Ho ^ Mi and Do = D\), then the firing rate changes continuously with time constant | = 1. 
If, instead, we change D and keep fi constant (by setting jio = Mi and Do ^ D\), the firing rate has jump discontinuities 
at time and fo, and changes with a faster time constant of = |. 



B.2.2 The time dependent cross-covariance 

We now look at the effects of changes in the input parameters on the conditional firing rate, Vy\w ( T ; f )• We first derive 
look at lag T = 0. The quantity Vy|v^(0,f) quantifies the tendency of the neurons to fire together. The conditional 
distribution, P(t,V \ W(t) — 1), of V(t) given that W crossed threshold at time t is a Gaussian with mean and variance 
given respectively by 

m c (0,t) = m{t)+p(t){\ —m(t)) 

and 

o 2 {Q 1 t) = o 2 {t){l-p(t)) 

where p(t) — is the sub-threshold correlation and m(t), <J 2 {t), and 7(f) are derived in the previous subsection. 
The firing rate at lag x = is then given by 

v v]w {0,t) = v(m c {Q,t),a 2 (Q,t),D t ). 

We now derive the conditional firing rate for times t > and lags x > 0. We break the derivation into three cases. 
The distribution of V(t + x) conditioned on a spike in W at time t is a one dimensional Gaussian. If t + x < to, the 
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mean and variance of this Gaussian are given by 

m c (x,t) — e~ T m c (0,t) + (1 -e~ T Wi 1 r „ . 

If t € [f ,fo], but t + x > to, the mean and variance are 

m c (x,t) = e -(( r + T )-'o)m c (fo-f,f) + (l-e-«' +T )- f o))Alo \ rn , 

a c 2 (T,f) = e - 2 «'+ T )-*)(T t 2 ( f o-f,0 + (l-^ 2(( ' +T) - ro) Po r G i°< f °j,f + T> f0 . 

Finally, when f > to, the mean and variance are 

m c {x,t) = e~ x m c (0,t) + (l-e~ T )no 

d 2 (T, f ) 

The conditional firing rate is then given by 



<J?(x,t) =e- 2 *oH0,t) + (l-e- 2 *)D Q ( t>t0 - 



V V \w(?,t)=J{m c (x,t) 7 (J 2 (x,t),D t+T ). 

Acknowledgements 

This work was supported by NSF Grants DMS-0604429 and DMS-0817649 and a Texas ARP/ATP award. 

References 

F. Apfaltrer, C. Ly and D. Tranchina, Population density methods for stochastic neurons with realistic synaptic kinetics: 
Firing rate dynamics and fast computational methods, Network -Comp Neural, 17, 373-418 (2006). 

L. Badel, W. Gerstner and M. Richardson, Transition-state theory for integrate-and-fire neurons, Computational and 
Systems Neuroscience 2010, Salt Lake City, UT. 

A. Bourgeat and M. Kern, Simulation of transport around a nuclear waste disposal site: the couplex test cases, 
Computational Geosciences (special issue), Springer (2004). 

C. H. Bruneau, F. Marpeau and M. Saad, Numerical simulation of the miscible displacement of radionuclides in a 
heterogeneous porous medium, Int J Numer Meth Fl, 49, 1053-1085 (2005). 

N. Brunei and P. E. Latham, Firing Rate of the Noisy Quadratic Integrate-and-Fire Neuron, Neural Comput, 15, 
2281-2306 (2003). 

Y. Burak, S. Lewallen, and H. Sompolinsky, Stimulus-dependent correlations in threshold-crossing spiking neurons, 
Neural Comput, 21 8, 2269-2308, (2009). 

A. N. Burkitt, A Review of the Integrate-and-fire Neuron Model: I. Homogeneous Synaptic Input, Biol Cybern, 95, 
1-19 (2006). 

R. Courant, K. Friedrichs and H. Lewy, Uber die partiellen Differenzengleichungen der mathematischen Physik, Math- 
ematische Annalen, 100, 32-74 (1928). 

P. Dayan and L. F. Abbott , Theoretical Neuroscience: Computational And Mathematical Modeling of Neural Systems 
Description, Cambridge, MA: MIT Press (2001). 

G. B. Ermentrout and N. Kopell, Parabolic bursting in an excitable system coupled with a slow oscillation, SIAM J 
Appl Math, 46, 233-253 (1986). 



20 



Galan, R. F., G. B. Ermentrout, and N. N. Urban. Stochastic dynamics of uncoupled neural oscillators: Fokker-Planck 
studies with the finite element method. Phys Rev E, 76(5), 561 10 (2007). 

C. W. Gardiner, Handbook of Stochastic Methods, Springer New York (1985). 

E. Godlewski and R A. Raviart, Hyperbolic systems of conservation laws In: Mathematiques et applications, Ellipses 
(1990). 

H. Hasegawa, Population rate codes carried by mean, fluctuation and synchrony of neuronal firings, Physica A, 388, 
499-513 (2009). 

L.M. Harrison, O. David and K.J. Friston. Stochastic models of neuronal dynamics, Phil Trans R Soc B, 360, 1075- 
1091 (2005). 

M. Iannelli, Mathematical theory of Age-Structured Population Dynamics, Mathematical Monographs, C.N.R. Pisa, 
(1994). 

J. Keener and J. Sneyd Mathematical Physiology, Springer Verlag, 2008 

P. Khorsand and F. Chance, Transient responses to rapid changes in mean and variance in spiking models, PLoS ONE, 
3, 1757 (2008). 

P. E. Kloeden and E. Platen, Numerical Solution of Stochastic Differential Equations, Springer New York, NY (1992). 

B.W. Knight, Dynamics of encoding in a population of neurons, J Gen Physiol, 56, 734-766 (1972). 

B. Lindner, Coherence and stochastic resonance in nonlinear dynamical systems, Humboldt University (2001). 

B. Lindner and L. Schimansky-Geier, Transmission of noise coded versus additive signals through a neronal ensemble, 
Phys Rev Lett, 86, 2934-2937 (2001). 

C. Ly and D. Tranchina, Spike train statistics and dynamics with synaptic input from any renewal process: a population 
density approach, Neural Comput, 21, 360-96 (2009). 

F. Marpeau, A. Barua and K. Josic, A finite volume method for stochastic integrate-and-fire models, J Comput Neu- 
rosci, 26, 445-57 (2009). 

N. Masuda, Simultaneous rate-synchrony codes in populations of spiking neurons, Neural Comput, 18, 45-59 (2006). 

I. Meda, I. Atwater, A. Bangham, L. Orci and E. Rojas, The topography of electrical synchrony among fi-cells in the 
mouse islet of Langerhans. Quart J Exp Physiol 69, 719D735, (1984). 

V. Melnikov, Schmitt trigger: A solvable model of stochastic resonance, Phys Rev E, 48, 2481-2489 (1993). 

D. Q. Nykamp and D. Tranchina, A population density approach that facilitates large-scale modeling of neural net- 
works: analysis and an application to orientation tuning, J Comp Neurosci, 8, 19-50 (2000). 

A. Omurtag, B. Knight, and L. Sirovich, On the simulation of large populations of neurons J Comp Neurosci, 8:1, 
51-63,(2000). 

S. Ostojic, N. Brunei, and V. Hakim, How Connectivity, Background Activity, and Synaptic Properties Shape the 
Cross-Correlation between Spike Trains, J Neurosci, 29, 33, (2009). 

P. Rasetarinera, Etude mathematique et numerique de la restauration biologique en milieux poreux, University Bor- 
deaux 1 (1995). 

A. Renart, N. Brunei and X. J. Wang, Mean-field theory of irregularly spiking neuronal populations and working 
memory in recurrent cortical networks In: Computational neuroscience: A comprehensive approach, ed by J. Feng, 
431-490(2007). 



21 



H. Risken, The Fokker-Planck equation: Methods of Solution and Applications, Springer- Verlag Berlin and Heidelberg 
GmbH & Co. K, (1989). 

J. de la Rocha, B. Doiron, E. Shea-Brown, K. Josic and A. Reyes, Correlation between neural spike trains increases 
with firing rate, Nature, 448, 802-806 (2007). 

E. T. Rolls, M. Loh, G. Deco, and G. Winterer Malone, Computational models of schizophrenia and dopamine modu- 
lation in the prefrontal cortex, PLoS Comput Biol, 9, 696-708 (2008). 

A. Renart, P. Song and X-J. Wang. Robust spatial working memory through homeostatic synaptic scaling in heteroge- 
neous cortical networks. Neuron, 38, 473 (2003). 

D. L. Ringach and B.J. Malone, The operating point of the cortex: neurons as large deviation detectors, J Neurosci, 
27, 7673-83 (2007). 

R. Rosenbaum, J. Trousdale, and K. Josic Pooling and correlated neural activity, Front Comput Neurosci 4:9 (2010). 

E. Salinas and T. Sejnowski, Impact of Correlated Synaptic Input on Output Firing Rate and Variability in Simple 
Neuronal Models , J Neurosci, 20(16):6193-6209, (2000). 

A. Schneider, T. Lewis, and J. Rinzel Effects of correlated input and electrical coupling on synchrony in fast-spiking 
cell networks. Neurocomputing 69, 1125-1129 (2006). 

M. Shadlen and W. Newsome, The Variable Discharge of Cortical Neurons: Implications for Connectivity, Computa- 
tion, and Information Coding , J Neurosci, 18:10, 3870-3896, (1998). 

E. Shea-Brown, K. Josic, J. de la Rocha and B. Doiron, Correlation and synchrony transfer in integrate-and-fire 
neurons: basic properties and consequences for coding, Phys Rev Lett 100, 108102 (2008). 

A. Sherman and J. Rinzel, Model for synchronization of pancreatic j3 -cells by gap junctions. Biophys J 59 547-559, 
(1991). 

G. Silberberg, M. Bethge, H. Markram, K. Pawelzik and M. Tsodyks, Dynamics of population rate codes in ensembles 
of neocortical neurons, J Neurophys, 91, 704-9 (2004). 

L. Sirovich, Populations of tightly coupled neurons: the RGC/LGN system, Neural Comput 20, 1179-210 (2008). 

H. Sompolinsky, H. Yoon, K. Kang and M. Shamir, Population coding in neuronal systems with correlated noise, Phys 
Rev E, 64(5), 051904(2001). 

T. Tchumatchenko, A. Malyshev, T. Geisel, M. Volgushev, and F. Wolf Correlations and synchrony in threshold neuron 
models, Phys Rev Lett 104(5), (2010). 

H.C. Tuckwell. Introduction to theoretic neurobiology, vol. 2. Cambridge University Press (1988). 

N. G. Van Kampen, Stochastic processes in physics and chemistry, North-Holland (2007). 

G. Webb, Theory of Nonlinear Age-Dependent Population Dynamics, Maecel Dekker, New York (1985). 



22 



